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This paper presents high fidelity Navier-Stokes simulations of the Boeing Subsonic Ultra Green Aircraft 
Research truss-braced wing wind-tunnel model and compares the results to linear MSC.Nastran flutter anal- 
ysis and preliminary data from a recent wind-tunnel test of that model at the NASA Langley Research Center 
Transonic Dynamics Tunnel. The simulated conditions under consideration are zero angle of attack, so that 
structural nonlinearity can be neglected. It is found that, for Mach numbers greater than 0.78, the linear flutter 
analysis predicts flutter onset dynamic pressure below the wind-tunnel test and that predicted by the Navier- 
Stokes analysis. Furthermore, the wind-tunnel test revealed that the majority of the high structural dynamics 
cases were wing limit cycle oscillation (LCO ) rather than flutter. Most Navier-Stokes simulated cases were also 
LCO rather than hard flutter. There is dip in the wind-tunnel test flutter/LCO onset in the Mach 0.76 - 0.80 
range. Conditions tested above that Mach number exhibited no aeroelastic instability at the dynamic pressures 
reached in the tunnel. The linear flutter analyses do not show a flutter/LCO dip. The Navier-Stokes simulations 
also do not reveal a dip; however, the flutter/LCO onset is at a significantly higher dynamic pressure at Mach 
0.90 than at lower Mach numbers. The Navier-Stokes simulations indicate a mild LCO onset at Mach 0.82, 
then a more rapidly growing instability at Mach 0.86 and 0.90. Finally, the modeling issues and their solution 
related to the use of a beam and pod finite element model to generate the Navier-Stokes structure mode shapes 
are discussed. 


I. Introduction 

The NASA Fixed Wing program has identified fuel efficiency and noise reduction targets as end products of 
funded research. To reach the goal of better fuel efficiency, propulsive and aerodynamic efficiency increases have been 
addressed in recent research efforts. The outcome of an N+3 phase 1 study conducted by the Boeing Company was 
the proposal to advance the Subsonic Ultra Green Aircraft Research (SUGAR) truss-braced wing (TBW) aircraft for 
further study. The advantages of the TBW configuration are its large aspect ratio, thin wings made possible by the 
added stiffness provided by a main and jury strut. The increase in aspect ratio results in a reduction of induced drag, 
while skin friction drag is reduced by maintaining laminar flow over thinner wings. One of the challenges of such a 
design is the likelihood of greater wing flexibility. This increased flexibility means that a multidisciplinary analysis of 
the vehicle aerodynamics, structure and control system is most likely required almost from the start of the design. The 
complexity of a wing/strut configuration designed to fly at a high subsonic Mach number suggests the use of higher 
fidelity tools to simulate the aerodynamic/structure interaction. 
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The Boeing Company is just now completing an N+3 phase 2 study of the SUGAR TBW concept. - This study 
included an analysis and design of a flight vehicle and the design and testing of an aeroelastic model in the NASA 
Langley Research Center Transonic Dynamics Tunnel (TDT). The Boeing designed SUGAR TBW model was tested 
in the TDT in December 2013 and January-February 2014. The purpose of that test was to assess the open and closed 
loop aeroelastic stability of the SUGAR TBW model in an unperturbed and oscillating free stream and to obtain a 
flutter boundary as a function of Mach number and dynamic pressure. Comparisons of the TDT wind-tunnel flutter 
onset with nonlinear MSC.Nastran flutter onset recently have been presented. The present paper will compare and 
contrast unsteady aeroelastic simulation results with the experimental data from the TDT test. For the present paper, 
the data of interest are the open loop flutter and limit cycle oscillation (LCO) onset boundaries in an unperturbed free 
stream. Low angle of attack conditions will be considered so that structural nonlinearity can be neglected. 

The primary focus will be simulations performed with the unstructured FUN3D Navier-Stokes computational fluid 
dynamics (CFD) code. Before unsteady aeroelastic results can be shown, the technique developed here to effect 
the fluid/structure data transfer will be discussed. The TDT wind-tunnel data to be presented in this paper will be 
discussed. The finite element model and flutter onset simulation results using the Boeing MSC.Nastran linear model 
of the SUGAR vehicle will be discussed. Finally CFD simulations of the vehicle aeroelastic response are presented 
and compared to the doublet lattice method (DLM) and experimental data. 


II. Wind-Tunnel Test Setup 


A. Transonic Dynamics Tunnel 

The Langley TDT is a unique national facility dedicated to identifying, understanding, and solving relevant aeroelastic 
and aeroservoelastic problems. The TDT is a closed circuit, continuous flow, variable pressure, wind-tunnel with a 
16 ft square test section with cropped corners. The tunnel uses either air or a heavy gas as the test medium and can 
operate at total pressures from near vacuum to atmospheric. It has a Mach number range from near zero to 1.2 and is 
capable of maximum unit Reynolds numbers of about 3 million per foot in air and 10 million per foot in heavy gas. The 
TDT currently uses 1,1, 1,2 tetrafluoroethane(R-134a) ’ as a heavy gas test medium. The TDT is specially configured 
for flutter testing, with excellent model visibility from the control room and a rapid tunnel shutdown capability for 
model safety. Testing in heavy gas has important advantages over testing in air: improved model to full-scale similitude 
(which results in heavier, easier-to-build models with lower elastic mode frequencies), higher Reynolds numbers, and 
reduced tunnel power requirements. Results presented are at nominal angles of attack in the range — 1 < a < 1, Mach 
numbers in the range 0.70 <M< 0.80 corresponding to Reynolds numbers of roughly 1 million to 2 million based on 
mean aerodynamic chord. 

B. Wind-Tunnel Model 

The SUGAR TBW wind-tunnel model is a 15 percent semispan side wall mounted model. The SUGAR TBW concept 
was developed by The Boeing Company and is described in reference 1. The wing and mount system for the SUGAR 
TBW test were designed and fabricated by NextGen Aeronautics in Torrance, California. NextGen Aeronautics was a 
subcontractor to The Boeing Company. Details can be found in reference 2. The model includes main and jury struts, 
pylon and annular nacelle, wing and fuselage. The vehicle empennage is not included in the wind-tunnel model. The 
fuselage includes an offset to bring the entire fuselage/wing outside the wind-tunnel wall boundary layer. The wing and 
struts are composed of spars for load bearing and pods to provide aerodynamic shape. The model is aeroelastically 
scaled for testing in heavy gas which, in addition to the advantages enumerated in the previous section, minimizes 
actuator bandwidth requirements. 

Figure 1 shows a photo of the SUGAR TBW model installed in the TDT. The scaled fuselage length was 18.8 
ft, but it was shortened to 13.4 ft to reduce weight and potential torsional dynamics with the turntable mount. The 
mean aerodynamic chord is 16.5 inches. The model wing span to the center line is 12.75 ft (about 75 percent of 
the tunnel width) with an additional 2.25 inches added to the half fuselage, the approximate wind-tunnel boundary 
layer displacement thickness. The wing, strut and jury use a spar-pod design where the scaled stiffnesses (El and 
GJ distributions) are designed into the flanged aluminum spar, and the aerodynamic shape is provided by discrete 
fairings, or pods, mounted to the spars. The pods were made of graphite epoxy and attached to the spars using brackets. 
The nacelle and pylon are not structurally scaled, but have the correct mass and inertia. The fuselage consists of an 
aluminum internal structure with fiberglass skin panels. Lead weights are attached to the wing spar in several locations 
so the wing has the proper scaled mass. 
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One of the the objectives of the SUGAR TBW test was to demonstrate flutter suppression. As a result, the wing 
has two trailing edge control surfaces, one near the wing tip and one just outboard of the engine pylon. Each aileron is 
driven by a vane-type hydraulic actuator with position measured by a rotary variable differential transducer (RVDT). 
Other instrumentation included 11 strain gauges and 21 accelerometers. The model was supported by a turntable on 
the east wall of the TDT as shown in Figure 1, and the model angle of attack was varied between -3 and +5 degrees. 

III. Computational Model Development 

A. Finite Element Structural Model 

Several MSC.Nastran finite element models (FEMs) of the SUGAR TBW configuration were generated over the 
course of this project. A full-scale, full-span version with highly detailed wing structure FEM was used in the de- 
sign/optimization studies described in reference . As part of the aeroelastic wind-tunnel model development process 
several semi-span, beam-rod FEMs were also developed. This FEM was refined as the details of the model structural 
design matured. The scaled mass and stiffness are reflected in version 19 of the beam-rod FEM shown in Figure 2a. 
This FEM includes updates to produce a good correlation with in-tunnel ground vibration test (GVT) measured modal 
frequencies and mode shapes. The structural damping has also been estimated from GVT to be around one percent. 
Note that a version 20 FEM with a better correlation with GVT data is now available, but is not used in this study. 

Version 19 of the beam-rod FEM has approximately 400 nodes with 200 elements including CBEAM, CONM2, 
RBAR , and RBE2 element types. This model consists of beams for the wing and strut spars, engine strut and control 
surface hinges. Discrete weights such as control surface actuators, spar/wing and engine spar/wing attachments are 
connected to the various beams by rigid body elements. Structural displacements are thus produced at a fairly complex 
arrangement of isolated node points. Accelerometer and strain gauge locations consistent with the physical model 
are included. Figure 2b shows the doublet lattice aerodynamic box layout used for the linear flutter analyses. The 
aerodynamic model was linked to the structural model using a combination of infinite plate and infinite beam splines. 
There are four major spline regions for the four major aero groups seen in Figure 2b: engine/nacelle, wing, strut, and 
fuselage. Also, each control surface has spline points associated with them as well. Additional modeling details can 
be found in reference 2 and a version of this FEM was used in the study described in reference 3. Both the linear and 
high fidelity nonlinear simulations discussed here use the normal structural modes from a vehicle finite element model 
that does not include any pre-stress. Modal frequencies for the first 16 modes are shown in Table 1. It can be seen 
that the first five modes include a variety of in-plane and out-of-plane wing bending and wing/nacelle torsion. Higher 
frequency modes also include strut bending and strut torsion. 


B. High Fidelity Nonlinear Computational Fluid Dynamics Model 

The high fidelity simulations are performed by the FUN3D CFD code. The Navier-Stokes code FUN3D (fully un- 
structured Navier-Stokes three-dimensional) is a finite-volume unstructured CFD code. Flow variables are stored at 
the vertices of the grid. The compressible Reynolds-averaged Navier-Stokes (RANS) flow equations and the Spalart- 
Allmaras turbulence model are solved on a tetrahedral grid. The Venkatakrishnan flux limiter is used and flux con- 
struction is performed by the low dissipation scheme. 

The SUGAR TBW surface mesh used in the FUN3D calculations is shown in Figure 3. A semi-infinite tetrahedral 
volume grid having 4.5 million nodes was developed with the VGR1D software. The wind-tunnel wall is treated 
as a symmetry plane. As in the wind-tunnel model, an offset has been added to account for the effect of the tunnel 
wall boundary layer. Upstream, downstream and side boundaries are placed at 10 wing span lengths from the model. 
Normal spacing of the first grid point off the vehicle surface is at y + ~ 1 everywhere. 

Steady state and subiterative solutions are accelerated to convergence by the use of local time stepping. Domain 
decomposition is used to enable distributed parallel computing. The solutions are computed on 288 — 800 processors 
depending on the platform. 

The aeroelastic simulations in this paper use mode shapes internal to the FUN3D code coupled with the generalized 
equations of structural dynamics. For a moving mesh, the conservation equations are written in the Arbitrary Lagrange 
Euler (ALE) formulation. The mesh deformation is accomplished by treating the mesh motion as analogous to a linear 
elasticity problem. The material properties vary based on distance to the nearest solid boundary. Elements near a solid 
boundary are made significantly stiffer by specifying the value of the Youngs modulus (E). 
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C. Force and Deflection Transfer Between Dissimilar Structure/CFD Models 

A key step in performing aeroelastic analysis is the transfer of structural deflections and fluid forces between the 
structure and the CFD surface nodes. Methods have been developed based on various interpolation and finite element 
shape functions. To name a few examples in use, there is the inverse isoparametric mapping, third order spline, 
nonuniform B-splines or those based on quadric surface basis functions. In the case of the FUN3D fluid/structure 
interaction, the Discrete Data Transfer Between Dissimilar Models (DDTBDM) software ' is often employed. That 
software uses a derivative of the inverse isoparametric mapping algorithm presented by Farhat, Lesoinne, and LeTal- 
lec. This method works well when the target and source nodes are in reasonably close proximity. 

As discussed earlier, the SUGAR TBW model fuselage, wing, trusses and nacelle pylon are modeled structurally 
by a collection of flexible and rigid beams with discrete masses. To accomodate the fluid/structure transfer for these 
various elements, the Boeing MSC.Nastran SUGAR model defines four infinite plate splines (IPS) and the beam 
nodes and aerodynamic boxes associated with each of these splines. Because the FEM structure and DLM aerodynamic 
models have similar coarse node spacing, it is straightforward to define spline regions that encompass all the structure 
and aerodynamic nodes. Furthermore, in the current DLM aerodynamic model, the wing and strut and nacelle boxes 
are not connected (see Figure 2). Different splines can be defined without discontinuities appearing between spline 
regions. It is more challenging when the structure FEM and CFD node spacings widely differ. Particularly, when the 
continuous CFD surface mesh has very tight grid spacing between several splined regions, any discontinuities at spline 
region intersections will be transmitted to the CFD surface. Such discontinuities in a CFD surface mesh are, of course, 
not tolerated. 

Methods of interfacing beam models with CFD surfaces that use an extrapolation of beam nodal displacements and 
rotations have been in existence for some time. • Using this approach, it is possible to extrapolate beam deflections 
using rigid body rotations and translations to CFD nodes that are far from the associated beam node. This approach 
has been very successful for applications such as an isolated rotor blade or wing where the transfer is between a simple 
beam and its associated CFD surface. - The DDTBDM software also has this capability to transfer between a fairly 
simple beam model and a CFD surface mesh. 

The SUGAR TBW model is a challenge because of the complex interconnection of multiply oriented beams. 
Because of this complexity, existing methods for transfer of beam deflections to CFD surface meshes have been found 
to be of limited use. This has motivated the development of an alternate approach. The approach used here starts by 
defining regions similar to the IPS spline regions defined in the MSC.Nastran model. Six regions rather than four are 
defined because the wing is divided into inboard and outboard regions. The six spline regions are shown in Figure 
4a. Also, rather than IPS splines, sixth order polynomial functions are defined in local optimally oriented coordinate 
systems. The user is able to choose terms up to sixth order, and to constrain end displacements and or slopes in the 
local coordinate system. This means, for instance, that if a swept wing region has no chordwise flexure, a simple wing 
twist or bending can be chosen that has constant slope in the swept chordwise direction. The mathematical description 
of the displacement in the local coordinate system is obtained as an over-constrained problem solved by singular value 
decomposition. Once the modal deflection has been calculated in the local coordinate system, it is transformed back 
into the global coordinate system. 

Two additional steps in the transfer process are added. First, an intermediate grid with mesh spacing that is coarser 
than the CFD mesh and finer than the beam model is used in the transfer process. The intermediate grid model is 
shown in Figure 4b. This intermediate internal solid tetrahedral mesh was created using the VGRID software from 
the original SUGAR TBW CAD model, by reversing the orientation of the surface mesh. This step transfers the beam 
deflections, to the intermediate grid with node points at or near the CFD surface. This allows the use of the DDTBDM 
inverse isoparametric mapping interpolation for the final transfer to the CFD surface mesh. 

The intermediate grid required smoothing of grid points that are at the boundary of or outside of spline regions. 
This configuration presents an additional challenge in that the wing has a supercritical profile shape. Roughly the 
aft five percent of the wing is exceptionally thin. This means that a small error in the displacement can result in the 
upper and lower trailing edge surfaces crossing over each other. In addition, some of the mode shapes have quite large 
curvature. This amplifies the problem of grid crossover. A Laplace smoother, such as that used in the DDTBDM 
software, results in unacceptable mode shape displacements in such areas of high curvature and thin body thickness. 
The solution was to go to a quadric surface based smoother. Rather than smoothing toward zero curvature using a local 
node support, a quadric surface smoother defined using a local cluster of grid points will retain a non-zero curvature. 
Furthermore, the smoothing is only targeted to small areas mainly between spline regions, omitting areas that require 
no smoothing. The result is a smooth deflected surface without slope discontinuities that also helps portions of the 
wing that are very thin, but also have very high curvature, retain the correct modal curvature. 

The smoothed intermediate grid deflections are transfered to the CFD surface mesh using the DDTBDM software. 
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It was also found that the inverse isoparametric interpolation would produce, in areas of high curvature (e.g. in some 
strut modes), slightly perceptible triangular faceting corresponding to the triangular surfaces of the intermediate grid 
tetrahedra. Once again, targeted smoothing was applied using the quadric surface smoother. Targeting small regions 
of the vehicle is intended to reduce as much as possible the degrading effect the smoothing has on the transfered the 
mode shape. 

The final product of the mode shape transfer can be seen in Figures 5- 10. Modes 3-5 are shown because these 
modes will be seen to be the main participants in the limit cycle oscillations (LCO) to be discussed. Mode 4 also 
has large nacelle pitching and translation. This presented challenges with smoothly connecting the nacelle and wing 
spline regions. Mode 9 is shown because it exhibits a high-curvature strut bending, as well as a very challenging 
connection between the deflected strut and the wing. A detail image of the strut/wing connection is shown in Figure 9 
to illustrate the smooth transition between spline regions in that area. Mode 15 is shown because of the complex wing 
bending and twist that it exhibits. Finally, the modal assurance criteria (MAC) was calculated using the mathematical 
representations defined in the six regions. The MAC was calculated by mapping the displacements back to the node 
locations of the original beam model. Table 1 shows the MAC for each mode up to 16. All the modes have MAC 
values at 92 percent or higher. The first 5 modes, that are mainly responsible for the LCO behavior to be discussed, 
have MAC values at 96 percent or better. 


IV. Results 


A. Linear Flutter Data 

For the linear flutter analysis, MSC.Nastran SOL 145 (PKNL method) was used. This solution procedure combines 
normal modes with a flutter analysis. For the PKNL method, a match set of densities and velocities are provided for 
each Mach number. For this analysis, the densities and velocities are calculated based on the R134a test medium. 
Figure 11 shows the V-g-f plot for Mach number 0.75. This is a typical plot seen for the range of Mach numbers 
analyzed for the wind-tunnel test. Here, it is noted that the flutter mechanism is a coalescence of modes 3 ( 2 nd out-of- 
plane) and 4 (1” torsion) at a dynamic pressure of approximately 59.5 psf and a frequency near 10 Hz. This analysis 
assumed zero structural damping; however, if 2 percent structural damping is assumed, the flutter crossing in terms of 
dynamic pressure is around 70 psf. The flutter onset for all Mach numbers is shown in Figure 12. As can be seen from 
that figure, the linear flutter analysis predicts onset at much lower dynamic pressures than the experimental flutter/LCO 
onset. 


B. Experimental Flutter and LCO Data 

For flutter testing, the goal is to establish the Mach number and dynamic pressure of flutter or limit cycle oscillation 
(LCO) onset. Operationally at the TDT, Mach number and dynamic pressure are not independent of each other and 
vary in unison with changes in the RPM of the TDT main drive motor at any given tunnel total pressure. A small 
portion of the heavy gas operating envelop of the TDT is shown in Figure 13a. Here, Mach number and dynamic 
pressure vary along lines parallel to the the total pressure (h) lines in psf from lower left to upper right as the main 
drive RPM increased. With the exception of several air checkout runs, all experimental data for the SUGAR TBW was 
acquired in heavy gas with a tunnel total pressure between approximately 250 and 450 psf. 

Typical flutter testing in the TDT relies on naturally occurring tunnel turbulence to perturb the model. Mach 
number and dynamic pressure are gradually increased along each h-line until a flutter point is obtained. Time history 
data are generally acquired along each h-line up to and including the flutter point. The determination of flutter is 
made visually by the test engineer, from digital strip chart observations, modal frequency coallescence or other visual 
assessments. Once a flutter point is obtained, Mach number and dynamic pressure are reduced to stabilize the model 
by gradually reducing main drive RPM or by opening the tunnel by-pass valves if warranted by the severity of the 
instability. Repeating this process over multiple h-lines allows the test engineers to define the flutter/LCO onset 
boundary. 

During the initial phases of testing, test engineers are more prone to trigger the by-pass valves at the first sign of an 
instability. As a test progresses, the test team becomes more familiar with both the aeroelastic response characteristics 
of the model and the effectiveness of the by-pass valves. For the SUGAR TBW model, points that were initially 
identified as flutter were later determined to be LCO or intermittent LCO. The experimental crossing appeared to be 
shallow and was perhaps sensitive to actual model damping. This can render the model amplitude and angle of attack 
dependent versus other wind-tunnel models that have steep crossings. 
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Most of the flutter testing for the SUGAR TBW test was conducted at angles of attack of —3, —1, 1, and 3 degrees. 
Flutter/LCO onset boundaries were found to vary with angle of attack. These data are presented in detail in reference 2. 
This angle of attack dependence was also noted in the analytical study described in reference 3. Because the present 
analysis focuses on lightly loaded conditions, so that structural nonlinearity can be neglected, only the minimally 
loaded experimental data at —1 and + 1 degrees are discussed here. 

The flutter/LCO data for the SUGAR TBW model at —1 and + 1 degrees are shown in Figure 13a and 13b, 
respectively. As with the MSC.Nastran flutter analysis, the experimental flutter mechanism was a coalescence of the 
2 nd out-of-plane bending and l If torsion modes at a frequency near 10 Hz. In post test data analysis, flutter/LCO 
onset was identified by a coalescence to a well defined single power spectral density (PSD) peak at around 10 Hz. 
The assessment of stability versus weak or strong LCO was made based on test observations, examination of the 
time histories (e.g. LCO amplitude) and PSDs of the acquired data. It should be noted that the SUGAR TBW wind- 
tunnel test was recently completed, so these assessments should be viewed as preliminary. The red symbols represent 
flutter or strong LCO data points, the blue symbols represent weak LCO points, and the open symbols represent stable 
points. Data points where flutter suppression control laws were opened to assess stability are included on this plot. 
This approach was helpful in establishing part of the backside onset boundary and for filling in flutter points within the 
unstable region. For + 1 degrees, the model exhibited some sustained torsional response at Mach numbers near 0.76. 

Tables 2 and 3 provide data for a few of the test points plotted in Figure 13. Test point numbers along with 
Mach number and dynamic pressure identify the conditions. The absolute maximum half peak-to-peak magnitude 
of Z-acceleration are tabulated for the forward wing tip and nacelle pylon accelerometers. The absolute maximum 
half peak-to-peak amplitude is used as a measure because, at many conditions, the LCO was intermittent, or it was 
necessary to open the by-pass valve to protect the model. In these cases, RMS values or average peak values would 
have been an inaccurate estimate of the largest LCO amplitudes. The reason for the accelerometer choice is that these 
accelerations are the best indicators of the participation of the modes primarily responsible for the LCO. The wing 
tip Z-acceleration provides an indication of the level of participation of the 2 nd out-of-plane mode. The nacelle Z- 
acceleration is an indication of the \ st torsion mode participation. These data will be used to make comparisons with 
FUN3D LCO amplitudes. There are some Mach number/dynamic pressure combinations that are repeated in these 
tables. It can be seen by looking at those repeat conditions that there is considerable scatter in the peak LCO amplitude 
observed. Additional analysis will most likely be required to fully understand the LCO amlitudes shown here. Finally, 
based on the LCO data points in Figure 13, a projected experimental flutter/strong LCO and weak LCO onset lines 
has been estimated and plotted in Figure 12 along with the flutter/LCO onset from the linear and the Navier-Stokes 
analyses. The projected onset lines are mean values of the — 1 and + 1 degree data. It is clear that there is a considerable 
difference in the dynamic pressure between what was considered weak and strong LCO onset especially in the lower 
Mach number range. By comparing these estimated LCO onset lines with the data plotted in Figure 13 it will be seen 
that some conditions within the weak LCO onset envelop were considered stable conditions. The estimated weak LCO 
onset line in the Mach 0.70 — 0.75 range, therefore, has considerable uncertainty associated with it. 

C. High Fidelity Nonlinear Computational Fluid Dynamics Flutter and LCO Data 

Unsteady aeroelastic simulations have been performed with the Navier-Stokes code FUN3D. The solutions are com- 
puted using the linear modes created by the process discussed in a previous section of this paper. The structural 
damping has been set to g - 0.02 (or damping ratio £ = 0.01). All of the conditions simulated here are at zero degrees 
angle of attack. Before discussing the flutter and LCO behaviors of the SUGAR TBW, a study of the time step size, 
CFL number, and number of iterations necessary for accurate solutions was performed. Figure 14 shows the devel- 
oping LCO for a condition at Mach 0.82. Three nondimensional time step sizes (4.0, 2.6, 1.7) were used to assess 
convergence. Each successive time step is a factor of 0.65 the previous time step size. The number of time steps per 
third mode cycle are shown in Table 4. Based on the LCO amplitudes seen in Figure 14, it is clear that the smallest 
time step is very near the converged LCO amplitude. The flow field CFL number used is 75, and the turbulence CFL 
number is 50, with 17 — 20 subiterations per time step. Figure 15 presents a random sampling of 100 lift and moment 
coefficient subiteration histories. The values plotted are actually ((;„ — ^i )/(^20 — ^i) for the n th subiteration with 
£ = {Cl, Cm) • The lift and moment coefficients are typically quite well converged in that range of subiterations. These 
results give confidence that the results presented here are well converged for the current grid. Note that at this time, a 
spatial grid convergence study has not been conducted. 

The computing procedure is that steady state rigid solutions are obtained, then static aeroelastic converged solu- 
tions are computed with a prescribed structural damping ratio set at £ = 0.999. The time-accurate aeroelastic solution, 
with the correct structural damping, is commenced by giving all 16 modes an initial modal velocity. The solutions 
were continued until it was clear that an absolutely unstable solution, or a steady state LCO amplitude, or zero LCO 
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amplitude was evolving. Many solutions required over 20-30 cycles of LCO oscillation, to reach a converged ampli- 
tude. 

Figure 1 2 presents a summary of the conditions that were simulated with FUN3D. The open symbols identify 
stable conditions (no LCO). Many of the closed symbols represent cases that reached an LCO. Only five conditions 
were approaching what appeared to be an absolute instability, namely the three conditions (with closed symbols) at 
Mach number 0.86 and two at Mach 0.90. Two important observations can be made from this figure. First, it is clear 
that FUN3D is producing predominantly LCO’s rather than flutter. The predominance of LCO is consistent with the 
behavior of the model in the wind-tunnel, which saw very few cases that were flutter. In many cases the wind tunnel 
produced moderate to large LCO amplitudes. Because the structural model in FUN3D is linear, the LCO is produced 
entirely by nonlinear aerodynamics. The second observation is that the LCO onset dynamic pressures in the FUN3D 
simulations for Mach numbers between 0.75 and 0.80 are much closer to the experimental onset than the MSC.Nastran 
flutter onset. The wind-tunnel data, with the model at a = ±1 degree, show a sharp rise in the onset dynamic pressure 
just above Mach 0.80. In that Mach range, the Navier-Stokes simulations, with the model at a - 0 degree, produce 
only a slight increase in the dynamic pressure at LCO onset, and a sharp increase in the onset dynamic pressure 
between Mach 0.86 and 0.90. The wind-tunnel test did not produce flutter above Mach number 0.80, mainly because 
of the difficulty in traversing the flutter/LCO dip at the higher dynamic pressures at which it would have occured. 

How close are the simulated LCO amplitudes to the experimental LCO amplitudes? As seen in Figure 12, in the 
Mach number range 0.70 — 0.74, the FUN3D LCO estimated onset line is significantly below even the wind tunnel 
test weak LCO range, but passes through the test weak LCO range at Mach 0.78. The maximum peak Z-accelerations 
from the selected wind-tunnel conditions in Tables 2 and 3 can be compared with the simulation data in Table 5. 
Consistent with the lower onset dynamic pressures in the Mach 0.70 — 0.75 range, the computed LCO amplitudes tend 
to be larger than the amplitudes observed in the wind-tunnel test. In the Mach number range 0.78 — 0.80, the computed 
amplitudes are similar to the experimental LCO amplitudes. This can be verified in Figure 13 and Tables 2, 3 and 5. 

The test point number 3443 in Table 3 is a condition that had a burst of LCO with amplitude very similar to that 
computed by FUN3D. The experimental Mach number is 0.786 and the dynamic pressure is 97.7 psf. The wind-tunnel 
model forward wing tip Z-acceleration experienced a large amplitude burst, as seen in the time history in Figure 16. 
At the time of the burst, the by-pass valve was opened to protect the model. This accounts for the rapid drop in LCO 
amplitude. One telling feature of the LCO in this figure is the change in frequency content before and after the by- 
pass valve was opened. Before the opening, the oscillation occurs at a well defined frequency at about 10 Hz. After 
the by-pass valve opening, and the amplitude drops considerably but there is a broad frequency spectrum of model 
oscillation. The broad spectrum oscillation and its amplitude is indicative of model response to tunnel turbulence. The 
time range 0—1 second, over which the highest LCO amplitude occurs in Figure 16, is compared with the computed 
time history of wing tip acceleration. That comparison is shown in Figure 17. The computational results were obtained 
at a Mach number of 0.78 and dynamic pressure of 90 psf. The computed LCO amplitudes compare very favorably 
with the experimental data in that figure. Note that what appears to be a phase shift between the two results is due to 
the different predicted frequency than that of the experimental data. The dominant LCO frequency in the simulated 
results is higher than that of the experiment. This can also be seen in the spectral transform data in Figure 18. This 
aspect of the simulations is currently not well understood. Additional simulations with finer grids and/or an updated 
structural model may be required to fully assess the accuracy of the simulations. 

The computed LCO onset above Mach 0.80 reveals an interesting feature that can be observed in Figure 12 and 
in Figures 19-21. After the slight rise in onset dynamic pressure through the Mach 0.75-0.80 range, there is a mild 
LCO onset at Mach 0.82. See Figure 19. By observing the time histories in Figure 20, it can be seen that while LCO 
onset at Mach 0.86 occurs at a slightly higher dynamic pressure, instability grows more rapidly at that Mach number. 
The LCO onset at Mach 0.90, between 100-115 psf, is much higher than at lower Mach numbers. Mach 0.90 and 130 
psf also reveals the explosive LCO amplitude growth seen at Mach 0.86. The simulations at Mach 0.86, 1 15 psf and 
Mach 0.90, 130 psf were terminated with negative volumes. Accelerations at the end of the simulations were 80-110 
g’s and growing. 


V. Concluding Remarks 

This paper has focused on the simulation of the flutter and limit cycle oscillation of the Boeing Subsonic Ultra 
Green Aircraft Research truss-braced wing wind-tunnel model. Modeling issues and their solution related to the 
transfer of the truss-braced wing wind-tunnel beam and pod finite element structural mode shapes to the CFD surface 
have been discussed. The SUGAR TBW model is a challenge because of the complex interconnection of multiply 
oriented beams. Existing methods for transfer of beam deflections to CFD surface meshes have been found to be of 
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limited use. An alternate method of transfering modal deflections from the finite element model to the CFD surface 
mesh has been developed. 

The computed flutter and LCO at model angle of attack of 0 degrees have been compared with data from a recent 
wind-tunnel test of the Boeing Subsonic Ultra Green Aircraft Research truss-braced wing wind-tunnel model at the 
NASA Langley Research Center Transonic Dynamics Tunnel. The wind-tunnel test and model are discussed. The 
wind tunnel model test data is for angle of attack of ±1 degrees. The high fidelity Navier-Stokes LCO simulations 
have also been compared to linear flutter onset. It is found that the linear analysis predicts flutter onset dynamic 
pressure well below the wind-tunnel test LCO onset. Furthermore, the majority of the wind-tunnel test high dynamics 
cases were wing and nacelle LCO. The Navier-Stokes simulations generally also exhibit LCO rather than hard flutter. 
Both the linear and Navier-Stokes simulations correctly predict that the LCO is produced by a coalescence of the third 
and fourth structural modes. There is a transonic dip and then sharp rise in the wind-tunnel test LCO onset in the Mach 
0.78-0.80 region. Wind-tunnel conditions above that Mach number exhibited no aeroelastic instability at the dynamic 
pressures reached in the tunnel. The linear analysis does not show this rise. The Navier-Stokes simulations, however, 
show a very slight rise in LCO onset from Mach 0.71 to 0.86 and a sharp rise from 0.86 to 0.90. At Mach 0.90, the 
simulations were damped at 100 psf and below. An interesting feature of the Navier-Stokes simulations is the mild 
LCO onset at Mach 0.82, but a more explosive LCO onset at Mach 0.86 and 0.90. 

There are many questions that remain unanswered from the wind-tunnel test and the computations. The physical 
mechanism for the appearance of torsional modes at Mach 0.76 and + 1 degree angle of attack and not at other 
conditions remains unclear. Additional work is required to determine the physical mechanism causing the difference 
in LCO onset seen in the computations in the Mach number range 0.82 — 0.90. The wind-tunnel data exhibit an 
angle of attack dependency, indeed, the shape of the transonic LCO onset envelop may be angle of attack dependent. 
Previously published flutter analyses using linear aerodynamics and nonlinear structural mode shapes show the similar 
angle of attack dependency as the wind-tunnel data. The present Navier-Stokes simulations have been computed at 
zero degree angle of attack. Future high fidelity simulations are needed to asses how LCO/flutter onset is influenced by 
angle of attack. The influence of structural nonlinearity, in particular in-plane strut loading on flutter/LCO onset has 
not been investigated with the high fidelity CFD code. Finally, an investigation of the influence of CFD grid resolution 
and a repeat of the present computations using the version 20 FEM, that may have an improved GVT correlation, is 
warranted. 
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Table 1. Modal data 


Mode Number 

Frequency, Hz 

Description 

Modal Assurance Criteria 

1 

5.13 

I s ' out-of-plane wing bending 

0.999 

2 

6.03 

I st in-plane wing bending 

0.966 

3 

9.17 

2 nd out-of-plane wing bending 

0.968 

4 

11.3 

l v ' wing/nacelle torsion 

0.988 

5 

18.5 

wing bending 

0.976 

6 

25.9 

wing/strut torsion 

0.989 

7 

27.7 

wing/nacelle torsion/bending 

0.981 

8 

28.2 

wing/strut torsion 

0.925 

9 

29.1 

strut bending 

0.964 

10 

38.6 

wing/strut torsion/bending 

0.963 

11 

40.5 

wing/nacelle torsion 

0.934 

12 

46.9 

wing/strut torsion/bending 

0.969 

13 

48.6 

wing torsion/strut bending 

0.961 

14 

61.6 

wing torsion/bending 

0.960 

15 

68.4 

wing torsion/bending 

0.951 

16 

76.4 

wing torsion 

0.952 


Table 2. Experimental Absolute Peak LCO amplitudes for the a = — 1 degree data points 


Point No 

Mach No 

Dynamic Press (psf) 

Forwd Wing Tip Z-Acc (g’s) 

Nacelle Pylon Z-Acc (g’s) 

4076 

0.699 

90.5 

5.5 

0.6 

3900 

0.719 

81.8 

4.8 

0.5 

3557 

0.720 

96.0 

9.5 

0.9 

3578 

0.740 

100.7 

11.1 

0.9 

4082 

0.741 

99.5 

9.5 

0.9 

4084 

0.742 

99.8 

16.0 

1.4 

3169 

0.759 

80.7 

11.1 

1.0 

3384 

0.760 

92.4 

13.1 

1.6 

3387 

0.761 

93.0 

13.9 

1.8 

3187 

0.779 

84.4 

11.6 

1.0 

3198 

0.782 

85.1 

13.3 

1.0 
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oint N 

2161 

2165 

2170 

2171 

2172 

2175 

2784 

3599 

3212 

3423 

3624 

1995 

2796 

3940 

4150 

3443 


Table 3. Experimental Absolute Peak LCO amplitudes for the a = + 1 degree data points 


Mach No Dynamic Press (psf) Forwd Wing Tip Z-Acc (g’s) Nacelle Pylon Z-Acc (g’s) 


0.708 

82.7 

8.9 

1.0 

0.718 

84.8 

7.4 

0.9 

0.739 

89.1 

9.4 

1.0 

0.743 

90.0 

11.4 

1.1 

0.744 

90.1 

10.1 

1.0 

0.758 

93.0 

11.2 

0.9 

0.761 

84.0 

4.6 

0.8 

0.761 

106.0 

8.5 

0.8 

0.779 

85.1 

15.3 

1.2 

0.779 

97.7 

12.6 

1.5 

0.779 

110.4 

22.0 

2.0 

0.781 

73.1 

13.3 

1.0 

0.780 

87.4 

9.7 

0.8 

0.781 

94.4 

11.9 

1.3 

0.781 

109.7 

6.3 

0.7 

0.786 

97.7 

30.7 

2.4 


Table 4. Time step size and number per 
third mode cycle 


Df 

T. S. per 3 rd mode cycle 

4 

170 

2.6 

277 

1.7 

423 
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Table 5. FUN3D computed LCO amplitudes, a = 0 degree 


Mach No 

Dynamic Press (psf) 

Forwd Wing Tip Z-Acc (g’s) 

Nacelle Pylon Z-Acc (g’s) 

0.71 

70 

12 

1.5 

0.71 

80 

26 

3.0 

0.71 

90 

40 

4.0 

0.78 

70 

5 

0.3 

0.78 

80 

13 

0.3 

0.78 

90 

30 

0.7 

0.82 

75 

7 

0.4 

0.82 

85 

10 

0.6 

0.82 

100 

20 

1.0 

0.82 

120 

35 

2.5 

0.86 

85 

20 

0.8 

0.86 

95 

55+ 

2+ 

0.86 

105 

70+ 

3 

0.86 

115 

80+ 

3 

0.90 

115 

50+ 

6+ 

0.90 

130 

110+ 

12+ 
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Figure 1. TDT SUGAR model 
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(a) Boeing MSC.Nastran FEM of wind-tunnel model 



(b) Aero boxes used in linear flutter analyses 


Figure 2. MSC.Nastran structural and aerodynamic models 



Figure 3. CFD surface mesh 
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(a) Regions for functional representation of displacement. 



(b) Intermediate tetrahedral model. 


Figure 4. Fluid/Structure transfer process 
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Figure 5. Mode 3 (Beam - deflected (red) - undeflected (black), CFD surface - deflected (blue) undeflected (gray)) 



Figure 6. Mode 4 (Beam - deflected (red) - undeflected (black), CFD surface - deflected (blue) undeflected (gray)) 



Figure 7. Mode 5 (Beam - deflected (red) - undeflected (black), CFD surface - deflected (blue) undeflected (gray) 
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Figure 8. Mode 9 (Beam - deflected (red) - undeflected (black), CFD surface - deflected (blue) undeflected (gray)) 



Figure 9. Mode 9 close up (Beam - deflected (red) - undeflected (black), CFD surface - deflected (blue) undeflected (gray)) 



Figure 10. Mode 15 (Beam - deflected (red) - undeflected (black), CFD surface - deflected (blue) undeflected (gray)) 
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Figure 11. MSC.Nastran SOL 145, PKNL method, V-g-f plot for Mach 0.75 in R134a 
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Figure 12. Experimental and simulated flutter/LCO conditions 


19 of 24 


American Institute of Aeronautics and Astronautics 


Dynamic pressure, psf Dynamic pressure, psf 



(a) a = -1.0. 



(b) a = +1.0. 


Figure 13. TDT experimental flutter/LCO points at two angles of attack 
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Figure 14. Time step convergence, Mach 0.82 




(a) Selection of Cl sub-iterative histories (b) Selection of Cm sub-iterative histories 

Figure 15. Typical sub-iterative convergence 
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Figure 16. LCO time history, forward wing tip z-acceleration, experiment - Point 3443 



Figure 17. Comparison of LCO time histories, forward wing tip z-acceleration, computed - Mach 0.78, q 90 psf, experiment - Point 3443 
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Figure 18. Comparison of LCO frequency spectra, forward wing tip z-acceleration, computed - Mach 0.78, q 90 psf, experiment - Point 
3443 



Figure 19. FUN3D computed LCO time histories, forward wing tip z-acceleration, Mach 0.82 
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Figure 20. FUN3D computed LCO time histories, forward wing tip z-acceleration, Mach 0.86 



Figure 21. FUN3D computed LCO time history, forward wing tip z-acceleration, Mach 0.90, Dynamic pressure 130 psf 
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